#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <cstdlib>
#include "txtpara.h"
#include "KDindex.h"
#include "directsampling.h"
#include "tplate.h"
#include <boost/random.hpp>
#include <ctime>
using namespace std;
double inline normalize(double y, double min, double max, double maxdiff)
{
    return (y-min)/(maxdiff)*255;
}
//Generating normally distributed random numbers
double sampleNormal(double mean, double sigma)
{
    //Establish Mersenne twister random number generator
    //Set seed using UNIX time
    static boost::mt19937 rng(static_cast <unsigned>(std::time(0)));
    //Select Gaussian random distribution
    boost::normal_distribution < double >norm_dist(mean, sigma);
    //Use the form of function to generate a random number generator
    boost::variate_generator < boost::mt19937&,
        boost::normal_distribution < double >>normal_sampler(rng, norm_dist);
    return normal_sampler();
}
int main(int argc, char* argv[])
{
    if (argc != 5)
    {
        cout << "Usage:\n" << argv[0] <<
            " param_file training_image_file template_file paramfile"
            << endl;
        return 0;
    }
    /* direct sampling parameters */
    int NR_neighbor;            //Maximum number of neighbors
    int NR_try;                    //Maximum number of search
    float tolerance;            //tolerance of difference
    int xsize, ysize, zsize;
    double percentage;
    int unitsize = 1;
    int maxdv = 20;
    char* config = argv[1];
    NR_neighbor = getparaINT("NR_neighbor", config);
    tolerance = getparaFLOAT("tolerance", config);
    xsize = getparaINT("xsize", config);
    ysize = getparaINT("ysize", config);
    zsize = getparaINT("zsize", config);
    maxdv = getparaINT("maxdv", config);
    percentage = getparaFLOAT("percent", config);
    std::vector < coord > sample;
    //Step 1.  Generate sample points for the first time
    //>>>>>1.1 calculating Number of points
    int tol_sample = xsize * ysize * zsize * percentage;
    int N_sample = xsize * ysize * zsize * percentage;
    int interval = round(sqrt(xsize * ysize / N_sample));
    if ((xsize / (interval)) * (ysize / (interval)) > N_sample)
        interval++;
    //>>>>>1.2 generate random path
    srand(clock());
    std::vector < int >ids;
    int zpos = 0;
    int s_no = 0;
    for (int xpos = 0; xpos < xsize; xpos += interval)
    {
        for (int ypos = 0; ypos < ysize; ypos += interval)
        {
            coord s;
            s.x = xpos;
            s.y = ypos;
            s.z = zpos;
            sample.push_back(s);
            ids.push_back(s_no);
            s_no++;
        }
    }
    N_sample = sample.size();
    for (int xpos = 0; xpos < xsize; xpos++)
    {
        for (int ypos = 0; ypos < ysize; ypos++) {
            if (((xpos % interval) == 0) && ((ypos % interval) == 0))
                continue;
            coord s;
            s.x = xpos;
            s.y = ypos;
            s.z = zpos;
            sample.push_back(s);
            ids.push_back(s_no);
            s_no++;
        }
    }
    for (int i = N_sample; i < tol_sample; i++)
    {
        int ranid = N_sample + rand() % (xsize * ysize - N_sample);
        int tmp = ids[i];
        ids[i] = ids[ranid];
        ids[ranid] = tmp;
        //Take the rest of the random samples from the rest
    }
    //Random order
    for (int s_no = 0; s_no < tol_sample; s_no++)
    {
        int ranid = rand() % N_sample;
        int tmp = ids[s_no];
        ids[s_no] = ids[ranid];
        ids[ranid] = tmp;
    }
    N_sample = tol_sample;
    //Step 2. Perform the simulation for first fixed percentage of cells
    Index < L2 < float >>* R = createIndex();
    float* ti;
    int tilength, tiwidth, NR_variable;
    int start_id = 0;
    int level = 0;
    load_ti_tiff(&ti, argv[2], tilength, tiwidth, NR_variable);
    float* target = new float[int(xsize * ysize * zsize * NR_variable)];
    float* target_ = new float[int(xsize * ysize * zsize * NR_variable)];
    for (int zpos = 0; zpos < zsize * NR_variable; zpos++)
    {
        for (int ypos = 0; ypos < ysize; ypos++)
        {
            for (int xpos = 0; xpos < xsize; xpos++)
            {
                target[zpos * ysize * xsize + xpos * ysize + ypos] = -1;
            }
        }
    }
    for (int zpos = 0; zpos < zsize * NR_variable; zpos++)
    {
        for (int ypos = 0; ypos < ysize; ypos++)
        {
            for (int xpos = 0; xpos < xsize; xpos++)
            {
                target_[zpos * ysize * xsize + xpos * ysize + ypos] = -1;
            }
        }
    }
    NR_try = int(getparaFLOAT("Fraction", config) * tilength * tiwidth);
    Simulation_cont(R, ti, NR_neighbor, NR_try, tolerance, sample, ids,
        target, N_sample, tilength, tiwidth, 1, NR_variable, unitsize, start_id, level, xsize, ysize, zsize);
         
    //Step 3 Interpolation using regression
    //1. Read in template file tp.par
    tplate tp;
    load_template_3D(&tp, argv[3]);
    data_event_cont* dv = new data_event_cont;
    dv->tpos = new int[tp.tpsize];
    dv->pos = new int[2 * tp.tpsize];
    dv->value = new float[3 * tp.tpsize];
    //2. Read in model file model.txt  
    ifstream param;
    param.open(argv[4]);
    int nr_param;
    int add_param = 3;
    int tpsize;
    param >> nr_param;
    param >> tpsize;
    int k1,k2,k3;
    float* theta = new float[3 * tpsize * (nr_param + add_param)];
    for (int i = 0; i < 3 * tpsize; i++)
    {
        for (int j = 0; j < nr_param; j++)
        {
            param >> theta[i * (nr_param + add_param) + j];
        }
        if(theta[i * (nr_param + add_param) + 1]>0)
        	k1=0.0;
        else
        	k1=255.0;
        if(theta[i * (nr_param + add_param) + 2]>0)
        	k2=0.0;
        else
        	k2=255.0;
        if(theta[i * (nr_param + add_param) + 3]>0)
        	k3=0.0;
        else
        	k3=255.0;
        theta[i * (nr_param + add_param) + nr_param] = theta[i * (nr_param + add_param)]
        + k1 * theta[i * (nr_param + add_param) + 1]
        + k2 * theta[i * (nr_param + add_param) + 2]
        + k3 * theta[i * (nr_param + add_param) + 3];
        if(theta[i * (nr_param + add_param) + 1]<0)
        	k1=0.0;
        else
        	k1=255.0;
        if(theta[i * (nr_param + add_param) + 2]<0)
        	k2=0.0;
        else
        	k2=255.0;
        if(theta[i * (nr_param + add_param) + 3]<0)
        	k3=0.0;
        else
        	k3=255.0;
        theta[i * (nr_param + add_param) + nr_param + 1] = theta[i * (nr_param + add_param)]
        + k1 * theta[i * (nr_param + add_param) + 1]
        + k2 * theta[i * (nr_param + add_param) + 2]
        + k3 * theta[i * (nr_param + add_param) + 3];
        theta[i * (nr_param + add_param) + nr_param + 2] = theta[i * (nr_param + add_param) + nr_param + 1] - theta[i * (nr_param + add_param) + nr_param];
    }
    param.close();

    for (int ypos = 0; ypos < ysize; ypos++)
    {
        for (int xpos = 0; xpos < xsize; xpos++)
        {
            if ((target[xpos * ysize + ypos] - (-1)) < 1e-6)
            {
                //Not simulated
                //Find the neighbors
                get_dv_3D(dv, &tp, target, xsize, ysize, xpos, ypos, maxdv);
                if (dv->dvsize > 0)    //It has neighbors
                {
                    {
                        //Regression results were calculated using each neighbor
                        //Average regression results                   
                        float y;
                        float y_sum = 0;
                        float y_avg;
                        float q;
                        float var_sum = 0;
                        float var_avg;
                        int k = 0;
                        /*
                         * Solve the sum of y values y_sum and the sum of variance var_sum in all conditional data y=ax+b
                         * */
                        for (int i = 0; i < dv->dvsize; i++)
                        {
                            if (
                                    (0.01 > theta[(dv->tpos[i]) * 3 * (nr_param + add_param) + 5]) &&
                                    (0.01 > theta[(dv->tpos[i]) * 3 * (nr_param + add_param) + 6]) &&
                                    (0.01 > theta[(dv->tpos[i]) * 3 * (nr_param + add_param) + 7])
                               )
                                //The pvalue is too large and has no statistical significance. The later ones do not need to be read in
                            {
                                y = theta[(dv->tpos[i]) * 3 * (nr_param + add_param)]
                                    + theta[(dv->tpos[i]) * 3 * (nr_param + add_param) + 1] * (dv->value[i * 3])
                                    + theta[(dv->tpos[i]) * 3 * (nr_param + add_param) + 2] * (dv->value[i * 3 + 1])
                                    + theta[(dv->tpos[i]) * 3 * (nr_param + add_param) + 3] * (dv->value[i * 3 + 2]); //Obtain the value of each regression y of conditional data   
                                //The number obtained here is equivalent to compressing the result, so you need to decompress the result to 0-255            
                                q = theta[(dv->tpos[i]) * 3 * (nr_param + add_param) + 4];  //The square of the standard deviation is assigned to q
                                q = q * q;
                                y_sum += y;    //Find the y sum of all conditional data
                                var_sum += q;    //Sum the variance of all conditional data
                                k++;    //Find the number of conditional data                                      
                            }
                        }
                        y_avg = y_sum / k;    //Obtain the mean value of conditional data y_avg
                        var_avg = sqrt(var_sum) / k;    //Get sigma                         
                        {
                            target_[xpos * ysize + ypos] = y_avg;    //target_ temporary storage simulation point
                        }
                    }
                    {
                        //Regression results were calculated using each neighbor
                        //Average regression results 
                        float y;
                        float y_sum = 0;
                        float y_avg;
                        float q;
                        float var_sum = 0;
                        float var_avg;
                        int k = 0;
                        /*
                         * Solve the sum of y values y_sum and the sum of variance var_sum in all conditional data y=ax+b
                         * */
                        for (int i = 0; i < dv->dvsize; i++)
                        {
                            if (
                                    (0.01 > theta[((dv->tpos[i]) * 3+1) * (nr_param + add_param) + 5]) &&
                                    (0.01 > theta[((dv->tpos[i]) * 3 +1)* (nr_param + add_param) + 6]) &&
                                    (0.01 > theta[((dv->tpos[i]) * 3 +1)* (nr_param + add_param) + 7])
                               )
                                //The pvalue is too large and has no statistical significance. The later ones do not need to be read in
                            {
                                y = theta[((dv->tpos[i]) * 3 + 1) * (nr_param + add_param)]
                                    + theta[((dv->tpos[i]) * 3 + 1) * (nr_param + add_param) + 1] * (dv->value[i * 3])
                                    + theta[((dv->tpos[i]) * 3 + 1) * (nr_param + add_param) + 2] * (dv->value[i * 3 + 1])
                                    + theta[((dv->tpos[i]) * 3 + 1) * (nr_param + add_param) + 3] * (dv->value[i * 3 + 2]); //Obtain the value of each regression y of conditional data     
                                //The number obtained here is equivalent to compressing the result, so you need to decompress the result to 0-255           
                                q = theta[((dv->tpos[i]) * 3 + 1) * (nr_param + add_param) + 4];//Standard deviation squared q    
                                q = q * q;
                                y_sum += y;    //Find the y sum of all conditional data
                                var_sum += q;    //Sum the variance of all conditional data
                                k++;    //Find the number of conditional data                                       
                            }
                        }
                        y_avg = y_sum / k;    //Obtain the mean value of conditional data y_avg
                        var_avg = sqrt(var_sum) / k;    //Get sigma                         
                        {
                            target_[ysize * xsize + xpos * ysize + ypos] = y_avg;    //target_ temporary storage simulation point 
                        }
                    }
                    {
                        //Regression results were calculated using each neighbor
                        //Average regression results 
                        float y;
                        float y_sum = 0;
                        float y_avg;
                        float q;
                        float var_sum = 0;
                        float var_avg;
                        int k = 0;
                        /*
                         * Solve the sum of y values y_sum and the sum of variance var_sum in all conditional data y=ax+b
                         * */
                        for (int i = 0; i < dv->dvsize; i++)
                        {
                            if (
                                    (0.01 > theta[((dv->tpos[i]) * 3 + 2)* (nr_param + add_param) + 5]) &&
                                    (0.01 > theta[((dv->tpos[i]) * 3 + 2)* (nr_param + add_param) + 6]) &&
                                    (0.01 > theta[((dv->tpos[i]) * 3 + 2)* (nr_param + add_param) + 7])
                               )
                            {
                                y = theta[((dv->tpos[i]) * 3 + 2) * (nr_param + add_param)]
                                    + theta[((dv->tpos[i]) * 3 + 2) * (nr_param + add_param) + 1] * (dv->value[i * 3])
                                    + theta[((dv->tpos[i]) * 3 + 2) * (nr_param + add_param) + 2] * (dv->value[i * 3 + 1])
                                    + theta[((dv->tpos[i]) * 3 + 2) * (nr_param + add_param) + 3] * (dv->value[i * 3 + 2]); //Obtain the value of each regression y of conditional data     
                                //The number obtained here is equivalent to compressing the result, so you need to decompress the result to 0-255        
                                q = theta[((dv->tpos[i]) * 3 + 2) * (nr_param + add_param) + 4];//Standard deviation squared q
                                q = q * q;
                                y_sum += y;    //Find the y sum of all conditional data
                                var_sum += q;    //Sum the variance of all conditional data
                                k++;    //Find the number of conditional data                                                                 
                            }
                        }
                        y_avg = y_sum / k;    //Obtain the mean value of conditional data y_avg
                        var_avg = sqrt(var_sum) / k;    //Get sigma                         
                        {
                            target_[2*ysize*xsize+xpos * ysize + ypos] = y_avg;    //target_ temporary storage simulation point 
                        }
                    }
                }
                else
                {
                    //The result is a random value between the maximum and minimum
                    target_[xpos * ysize + ypos] = target_[ysize * xsize + xpos * ysize + ypos]
                        = target_[2 * ysize * xsize + xpos * ysize + ypos] = (float)(rand() % 255);
                }
            }
        }
    }

    for (int zpos = 0; zpos < zsize * NR_variable; zpos++)
    {
        for (int ypos = 0; ypos < ysize; ypos++)
        {
            for (int xpos = 0; xpos < xsize; xpos++)
            {
                if ((target[zpos * ysize * xsize + xpos * ysize + ypos] - (-1)) < 1e-6)
                {
                    target[zpos * ysize * xsize + xpos * ysize + ypos] = target_[zpos * ysize * xsize + xpos * ysize + ypos];
                }
            }
        }
    }
    //Step 4 Output simulation results
    for (int zpos = 0; zpos < zsize * NR_variable; zpos++)
    {
        for (int xpos = 0; xpos < xsize; xpos++)
        {
            for (int ypos = 0; ypos < ysize; ypos++)
                cout << target[zpos * ysize * xsize + xpos * ysize + ypos] << " ";
            cout<<endl;
        }
    }
    return 0;
}
